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Abstract 

The analytic form of a new class of factorized Runge-Kutta-Chebyshev (frkc) stability polynomials 
of arbitrary order N is presented. Roots of frkc stability polynomials of degree L = MN are used to 
construct explicit schemes comprising L forward Euler stages with internal stability ensured through 
a sequencing algorithm which limits the internal amplification factors to ~ L 2 . The associated 
stability domain scales as M 2 along the real axis. Marginally stable real-valued points on the 
interior of the stability domain are removed via a prescribed damping procedure. 

By construction, frkc schemes meet all linear order conditions; for nonlinear problems at orders 
above 2, complex splitting or Butcher series composition methods are required. Linear order condi¬ 
tions of the frkc stability polynomials are verified at orders 2, 4, and 6 in numerical experiments. 
Comparative studies with existing methods show the second-order unsplit FRKC2 scheme and higher 
order (4 and 6) split frkCs schemes are efficient for large moderately stiff problems. 
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1. Introduction 


Runge-Kutta-Chebyshev methods are explicit numerical integration schemes with extended sta¬ 
bility domains derived from the optimality properties of Chebyshev polynomials 00- These 
methods are commonly applied to moderately stiff systems of semi-discrete equations of the form 

w' = f(t, w), (1) 


yielding an approximate solution w n at time t n = nT, defined on a spatial mesh of spacing h at 
points Xk, with Xk+i = Xk + h. Such systems arise naturally through application of the method 
of lines to parabolic systems. Runge-Kutta-Chebyshev methods may be broadly categorized as 
factorized or recursive in nature. 

Factorized Runge-Kutta-Chebyshev methods are formed from a sequence of forward Euler 
stages. These methods were first suggested by Saulev 0, Guillou and Lago jl] and were sub¬ 
sequently considered by Gentzsch and Schluter Q and van der Houwen 0. They have been 
applied at first-order and extended to second-order via Richardson extrapolation by various au¬ 
thors mum. Based on a strategy proposed by Lebedev [llj| , the dumka stability polynomials 
exist at orders 2, 3, and 4 121. 
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Recursive Runge-Kutta-Chebyshev methods were first described by van Der Houwen and Som- 
meijer [3] and rely on (three-term) recursion to generate a solution. They were introduced at 
second-order by van Der Houwen and Sommeiier 13] and, subsequently, other second, third, and 
fourth-order methods have been developed 13, HI, 0, E3, EH] • We note that alternative approaches 
with second-order accuracy involving Legendre polynomials have recently been proposed by Meyer 
et al. mm. At orders above 2 , for both factorized and recursive methods, composition tech¬ 
niques relying on Butcher series theory [Hj, [22| are typically used to satisfy the full set of order 


conditions [ 12 l 117 1 . 


This paper is organized as follows. In Section [2] the analytic form of the class of frkc sta¬ 
bility polynomials is presented. The construction of stable time-marching schemes based on the 
roots of these polynomials is outlined. Section [3] is given over to the derivation of the polynomial 
through consideration of associated recurrence relations. In Section Q] numerical tests are pre¬ 
sented confirming the order and efficiency properties of frkc methods. Conclusions are presented 
in Section [5] 


2. High-order factorized Runge-Kutta-Chebyshev 

2.1. General prescription 

Eq. Q] may be written in autonomous form by appending t to the vector of dependent variables 
for the system 

w ' = f(w). ( 2 ) 

Parentheses may be used in the remainder of this work to differentiate exponents from indices. We 
proceed by considering order N extended stability explicit Runge-Kutta schemes over L = MN 
stages 


L 

W L = W°+ Tj2aif(W l ~ 1 ), (3) 

i=i 

where W° = w n corresponds to the approximate solution w n at time level n, and W L yields «; n+1 
at a time T later. The timestep related to each stage is then given by r; = a{T. 

The frkc polynomial of rank TV, and degree L, is given by 


N 

Bti(z)=d» +2^df C kM (z), 

k =1 


( 4 ) 


where CkM denotes the the Chebyshev polynomial of the first kind of degree kM. The corresponding 
optimal real stability range is [—(3m, 0 ], where (3m = 2M 2 ctM, ctM = (jmN + 2)/3, and 72 « 0.87 
(with 7 m rapidly converging to 1 as M increases). In this limit, the polynomials generate 81%, 
74% and 73% of the optimal intervals for order 2, 4, 6 respectively (see Van Der Houwen [23] and 
Abdulle [H] for estimates of the optimal values for a*). The limiting step size is T = /3M/|A| m ax, 
where A are the negative-definite eigenvalues for the Jacobian of Eq. [2] We note that the form 
of Eq. 0] is consistent with the known result that Chebyshev expansions of stability polynomials 


27 1 


exist to arbitrary order [251 ]. Furthermore, following from a proposition by Lomax [261 ]. Riha 
confirmed the existence and uniqueness of optimal stability polynomials with L — N local maxima 
with value unity. A full derivation of the frkc polynomial expression given by Eq. [4] is provided in 
Section [3] 


2 










-1 


-0.5 


0 


-1 


-0.5 


0 


x //3 20 


x/fco 


Figure 1: Absolute values along the real axis for frkc stability polynomials corresponding to 
M = 20 at various values of N. Damped polynomials with uq = 0.05 are shown with solid lines in 
the left side panels, and for y > 0 in the right side panels; associated undamped polynomials are 
also shown with dashed lines on the left, and for y < 0 on the right. For N = 2: 720 = 0.9988, 
/? 20 = 1066.0; for N = 4: 720 = 1.02 1 5, /3 20 = 1623.9; for N = 6 : 7 20 = 1.0276, foo = 2177.5. 
Left side panels show plots of \R^j(x)\ ( x £ R): (a) |-R| 0 ( a; )l; (b) l-^ 2 o( a ')li ( c ) l^ 2 o( a ')l- Dotted 
lines indicate guide values at 1.0, 0.95, 0.0. Right side panels show \Rm\ = 1: (d) |i?| 0 | = 1; (e) 

l^ol = i; (f) l^2ol = i- 
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The order coefficients , which we refer to collectively as the order pattern , are determined 
by requiring that the (undamped) stability polynomial R%(z) = B^( 1 + z/M 2 cxm ), consisting of 
shifted Chebyshev polynomials, satisfies the linear order conditions 

RM in \ 0) = 1, n = 1, , TV. (5) 

This requirement is met by solving the TV-dimensional linear systenQ 


' ^>(1) 
- ^ (1) 


coupled with the constraint 


C 


(1) fl) 

nm\ l > 


M N ) (-1 'i 
° NM\ l ) -I L 


N 




d N 

a N 


( M 2 cxm ) 1 

(. M 2 a M ) N 


N 

< = i-W- 


( 6 ) 


(7) 


Following identification of the roots of the frkc polynomial the damped order TV scheme 

corresponding to Eq. [3] is determined via 


1 1 
M 2 a M 1 — cr 


( 8 ) 


In order to ensure a stable scheme for small perturbations from the real axis in the spectrum of 
Eq. [H it is necessary to introduce a suitable damping procedure. We find an effective prescription 
for the damped order TV scheme is given by 


_1_ 1 - m 

(1 — v)M 2 a M 1 — (1 — 2 W )Ci’ 


(9) 


where the damping is parameterized by the small positive quantity v, resulting in the real extent 
of the stability interval being reduced to (1 — v)[$m- The value of u = v 0 /N is regulated by means 
of the reference damping parameter t'o, such that maxima in \R\ along the real axis are scaled by 
approximately 1 — vq. 

For the case = 0.05, with M = 20, and for various values of TV, Fig. |T] illustrates the effect of 
the damping procedure. It is clear that the undamped polynomials are marginally stable at M — 1 
points on the interior of the stability domain along the real axis. (In fact, for sub-optimal olm, 
internal marginally stable points occur at M/2 — 1 locations for even values of M, or (M — l)/2 
locations for odd values of M.) Examples of the order patterns for M = 20 with TV = 2, 4, 6 are 
given in | Appendix A| 

The L-tuple [/q] has cardinality TV and regulates the implementation of damping in the scheme 
while preserving the nominal order of accuracy. The values of pi are obtained by tuning the damped 
stability polynomial Rm( z ) = n^ =1 (l + aiz) to meet the linear order conditions given in Eq. [5] We 
describe the procedure for the determination of the damping coefficients pi in Section 12.21 


lr The identity f (1) = rii=o((^^) 2 — * 2 )/(2* + 1) is useful here. 
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2.2. Identification of damping parameter L-tuple 

The elementary symmetric polynomial, cr™ = <j l <m IIi=i C?», is defined as the sum 

of all possible products formed from l unrepeated elements drawn from the first m elements of an 
L-tuple [ft]. The definition is extended by setting er™ = 1 and cr™ >m = 0. We associate the L roots 
Ci, in order of increasing real component 5R(0), with the damping coefficients /q by cycling through 
the N damping coefficients a total of M times. Newton-Raphson iterations then converge rapidly 
to the linear order conditions given by Eq. [5] The effects of the damping procedure on the stability 
domain are shown in Fig.Q] The stage intervals 77 given by Eq. [9] are complex in general, however, 
with dj = 0 , d\ = 1 / 2 , the standard first-order super-timestepping scheme [7, .§] is recovered with 
B\i = Cm- For TV > 1, either one or two values of t; have negative real parts. 

The presented prescription implements conjugate pairs separately thereby necessitating full 
complex arithmetic. Other than some penalty in the additional computational demand required, 
we find no practical disadvantage to preserving this model of treating each factor as distinct. We 
note that Lebedev Emu proposed a scheme which treats stages in pairs and, when applied to 
conjugate pairs, removes the need for complex arithmetic. 


2.3. Internal stability 

Schemes comprising a high number of stages are internally unstable if the sequencing of the 

1—1 1—1 n i3 HU- Lebedev 




stages is allowed to admit uncontrolled growth of numerical errors 
and Finogenov ,33] first suggested sequencing stages to manage uncontrolled growth of internal 
instabilities (see also (34J). Here, we present a straightforward algorithm for sequencing stages 
which limits the maximum amplification factor of internal instabilities to ~ L 2 . 

We define Vj^k = |1 + ajX k |, where Xk £ [—/3m, 0] are discrete values spanning the spectrum of 
Eq. |T2] The L-tuple [a;] is then ordered by holding the Li normed quantity 


( 1 

L \ 

max v hk , 

n v i k 

V=1 

3 = 1+1 ) 


1 


( 10 ) 


to a minimum value while l is increased from 1 to L. This procedure suppresses the growth of the 
internal stability functions Qj : k (x ) = II Lj I 1 +aix\, for j, k = 1, • • • , L, over x £ [—/3m, 0], and 
provides excellent internal stability properties with high numbers of stages at low computational 
cost. In Fig. [2J we plot the maximum internal stability function Q(x) = ma x.j } k(Qj,k(x)) for the 
test cases L ss 4000, 400, 40, with N = 2, 4, 6, and vq = 0.05. The optimization may be enhanced 
by concentrating the points Xk towards the bounds of the interval. (In this work a logistic function 
over a range of 15 is employed to generate the L sample points.) We observe the maximum internal 
amplification factor scales approximately as L 2 . Hence, the internal stability properties are well 
within the acceptable limits of modern computing precision for any practical problem. 

Consistent with these findings, we note that internal amplification factors of ~ 10 6 are quoted 
in the literature for rkc methods with 1000 stages 35|, and furthermore, a quadratic dependence 
on stage number is suggested by Sonnneijer et al. [151] - ROCK2 methods are reported to demonstrate 
amplification factors of ~ 10 9 at 200 stages by Hundsdorfer and Verwer [351] , suggesting internal 
instability growth rates 150 time larger than for rkc and FRKC2 schemes. 

We note that the Serk scheme is also limited in stage number, requiring 600 digits of precision 
for 320 stages, albeit principally due to severely ill-conditioned matrix systems used in calculating 
the stability polynomials by means of the Remez algorithm [T3]. A subsequent revision of the serk 
methodology has demonstrated a stability range which is four times larger (36| . 
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Figure 2: Maximum internal stability function Q(x) for L ss 4000 (upper lines), L sa 400 (middle 
lines), L sa 40 (lower lines): (a) N = 2 with M = 2000, 200, 20; (b) N = 4 with M = 1000, 100, 10; 
(c) N = 6 with M = 667, 67, 7. In all cases the default value of = 0.05 is used. Guidelines show 
values of L 2 . 
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3. Factorized Runge-Kutta-Chebyshev polynomial derivation 

We proceed by considering the one-dimensional diffusion equation 


dw d 2 w 
~dt = 


(11) 


The semi-discrete form of Ea. fTTl mav be written w' = h~ 2 Dw, where I? is a tridiagonal matrix with 
diagonal entries -2, subdiagonal entries 1, corresponding to a second-order central discretization 
of the spatial derivative. The eigenvalues of D are negative with a maximum magnitude of 4. 
Application of the numerical scheme given by Eq. [3] yields 


w 


n+1 


L 



W 


( 12 ) 


The frkc polynomial B^ may be derived by consideration of the canonical scheme given 
by Eci. [T2l over an extended timestep T, spanning time levels t n to t n+1 , and consisting of M seg¬ 
ments , with each segment comprising N stages. We write the solution state corresponding to w n as 
W°, and assume Wq = 1 and 0 = 0; more complex states may be constructed by superposition. 

The solution state corresponding to w n+1 is then obtained from W M = n Zi(I + h- 2 T l M D)W°. 
To aid the following discussion, Fig. [3] is provided to graphically represent solution states W m at 
different segment levels for the particular case N = 2. A reference point value of the solution state 
W M , at spatial index j. is shown as a black node. 

To proceed, we assume schemes consisting of m segments (m = 1, • • • , M — 1) are known which 
generate the solution states, W m = (I + h~ 2 Tj n D)W°. For m = 1, the solution state W 1 
spans 2N + 1 nodes from a given point profile W°. Successive states regenerate this pattern, but 
spanning ( 2mN + 1) nodes, with non-zero values interspersed by ( m/N — 1) zero valued nodes. 
We refer to the sequence of patterns over increasing values of m as a pattern flow. Illustrations of 
sample pattern flows are given in Fig. [3] 

Using Eq. [51 the components of the states W m may be recast as Wj' = Wj n ^”'7(1 ~ Q n ) ■ 
Over a single timestep, Eq. [12] then takes the simplified form 


mN 

w m = UDTw°, (is) 

i=i 

where D™ is a tridiagonal matrix with diagonal entries —([ n and subdiagonal entries 1/2. In terms 
of the elementary symmetric polynomials we have 



1=0 


(14) 


where c™ l are coefficients dependent on the scheme Eq. [T5] By induction, these coefficients have 
the properties 


L 0, m 


Ca 


j, l^m 


=(-i r 

4(< 


' 3 - 1,1 


+ c m - 

+ C j+1 


!■) 


r m —0 

c j^0, m u ’ 


C j, 2^0 


= — c 


' 3 , 1 - 1 * 


(15) 
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Figure 3: Graphical representation of construction of primitive recurrence relation between state value W j 
and solution states for m < L at intervals of N = 2; non-zero coefficients used in Eq. [18] corresponding to 
Gaussian polynomials [■£] are shown (P = 2 N — 1). Nodes at values of m which are non-integral multiples 
of N (empty circles) do not appear in the relation construction. Pattern flows emerging from three sample 
source distributions up to segment level m = 2 are shown (labeled a , 6, and c). Rays terminating from the 
filled node at m = L and originating at the apices of the sample distributions are shown summing to unity 
(solid lines). Rays which do not similarly project the solution pattern from m = 2 through to m = L sum to 
zero (dashed line). Also shown are the coefficients r^ k prescribed by Eq. [2D] at given values of k. 













The (to + l)-tuples [cq m ] fully determine W™ through the roots Q 1 of the associated polynomial 
B// defined by 

jdN 171 

< 16 > 

' ’ 1=0 

where (2) miv_1 is a normalization factor. Hence, the (TOfV)-tuple [r ; m ] is completely specified by 


3.1. Derivation 

The (to + l)-tuple [0 bi , er m ~ b , 0 bH ] is constructed from the elementary symmetric polynomials 
corresponding to the solution W™ , where a m ~ b indicates the ordered elements o-™~ b , ■ ■ ■, cr™“ b , 
zero superscripts denote multiplicity, and b = bv + bn. Through Eq. 1161 the (to + l)-tuple maps 
to the degree m - b L polynomial (- l) bL (z) bR (2)- ( ~ m - b '> N+1 B^_ b . Inserting [0 bL , a m ~ b , 0 bfl ] into 
Eq. [14] and appealing to the properties of the coefficients cj l l , as given by Eq. [TS] yields a direct 

correspondence to (—l) bi (^) bR Y / b g =o Cg)^T-bR+2g- Hence, we derive the association 



W 


■m—b 

j-bR+2g 


jdN 

(r,~\bR m-b 

' ( 2)( m ^ b ) JV “ 1 ’ 


(17) 


By Eq. [13l the solution state W m+1 generates a pattern scaled by a factor of 1/2 with respect to 
the pattern corresponding to the solution state W . Hence, a recurrence relation generating the 
correct pattern comprising any weighted average of (2) L ~ m W over available values of to will yield 
a valid solution state W . 

We define a ray as any connection on a uniformly spaced graph which passes through nodes 
on every segment level to, m = 1, • • ■ , M — 1. The sum of the recurrence weightings over any ray 
terminating at to = L must be unity if the ray originates at the origin of a pattern flow at to = 0, 
and zero otherwise. The coefficients of the Gaussian polynomials ? (fc = 1, • • • , P), denoted 

[■£] , possess the required properties. In Fig. [3] rays are shown summing to unity and zero, with a 
list of weightings satisfying these properties for all possible rays for the particular case of N = 2. 
Defining P = 2N + 1, the primitive form of the recurrence relation for W is 


N 


^=EH) W E 


G k rp y 


k= 1 


1=0 


i \ kN 

1 \ TTrL-kN 

W j-±G k +l 


(P-k)N 


—L-(P- k )N 
W 3-hGk+l 


i \ PN 

n wr N . 


(18) 

where G k = kP — k 2 is the degree of [^] for k < N. We note that the Gaussian polyno¬ 
mial [■£] possesses a unique representation as a summation of the binomial powers (1 + q 2 ) 9 , for 
9 = 0, ■ ■ ■ , G k / 2 , given by 

\G k 

^2r^ k qi Gk ~ 9 {l + q 2 ) 9 , (19) 


P 


L J 1 3=0 

where the coefficients r// k follow the generating function 


^^(-l) fc (2)^ fc (t) fe (^ 

k —0 g —0 


N 

0 l-t)l[(l + (t) 2 -2tc k ). 

k= 1 


( 20 ) 
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Then, using Eq. [19j we may recast Eq. [TS] in the form 


N g 

5=0 


fc= 1 


2=0 


kN 


„ r L-kN 
W j-g+2l 


( P-k)N 


w 


■L-(P-k)N 

3-9+21 


+ .i 


PN 


Wj 

( 21 ) 

Applying the association given in Eq. [T7] to the terms in Eq. [21] the recurrence relation for 
follows as 

N |Gfc 

-®M = E(~^ + ~ ^M-P+k] T r a’ k (^ z ) 9 + Bm-P- (22) 

/c—1 5=0 

We continue by noting that the generating function for BJ? derived from the recurrence relation 
given by Eq. [22] is 


E(*) fc < = 


EZoi^k 


k—0 


1 - E"=i (-!) fc+1 [( t) k - El=o rp k (2z)9 - (t)P 


(23) 


where bjf are coefficients determined by the seed states of . Appealing to Eq. [20 the generating 
function derived from the recurrence relation given by Eq. [50 is 


f> W = ZL + 2J2 , 5^2 ftr (1 - *) II( 1 + W 2 - 2tC ^ (24) 

k=o *=i 1 + W -2fb fe fc=1 

where 6^ are coefficients determined by the seed states of B^[. The normalization B£ (1) = 
bff (1) + J2k =l 2 b^ (1) has been imposed in order to fix the forms of the numerators in the separated 
fractions. 

Noting that the generating function for Ckm is = (1 — zt)/{ 1 + ( t ) 2 — 2 Ck), 

we conclude that Bj? = bff + ZYlZi b^Ckm- Consideration of the particular case N = M = 1 
indicates a correspondence between b\ and d\ is required in order to match the required solution 
pattern and normalization properties. A general correspondence between b^ and dff is established 
by considering successive values of N , with M = 1, for the limiting case dj? = 0, 0 < k < N. 
This completes the derivation of the analytic expression for the frkc stability polynomial given by 
Eq.H 


N 


4. Tests 

In this section numerical studies of two-dimensional two-species Brusselator diffusion-reaction 
problems are presented which confirm that high-order frkc stability polynomials meet all relevant 
linear order conditions and that the derived factorized numerical schemes are both stable and effi¬ 
cient. Split schemes, denoted frkCs, are obtained by means of complex splitting techniques: linear 
diffusion operators are treated via frkc methods, while nonlinear reaction terms are integrated us¬ 
ing standard Runge-Kutta techniques. The performance of the second-order accurate unsplit FRKC2 
scheme is compared to second-order rkc and CVODE2 codes . Finally, comparisons are presented 
of higher order split frkCs schemes (at orders 4 and 6), with fourth-order ROCK4 , and fifth-order 
CVODE schemes. 


PN 
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4-1. High-order splitting 

frkc stability polynomials satisfy linear order conditions to an arbitrary order of accuracy. This 
property may be exploited in solving numeric al p roblems for semi-linear stiff systems of equations 


through operator splitting methods [371 l38|, [39|, |40|. We note that in the literature, the linear 


and nonlinear terms of reaction-diffusion models have been decoupled under a variety of numerical 
integration techniques including: splitting methods [4l| . and previous references], Implicit-Explicit 
Runge-Kutta-Chebyshev (IMEX RKC) methods [ijE^], PIROCK :44|, and Local Linearization 
Runge-Kutta (LLRK) methods (45), 0] . 

High-order splitting has been shown to give rise to an order reduction effect in some reaction- 
diffusion cases [471 ]. For Dirichlet and Neumann boundary conditions, splitting techniques may 
result in order reduction at boundaries 35, 38]. It has also been observed that the full order is 


recovered on the interior of the computational domain when it is taken sufficiently far from the 
influence of the boundaries 48,0]. Boundary conditions for the separate operator updates are 


necessary to avoid order reduction effectively, however, as yet, no consistent treatment exists [50(. 

Assuming Eq. [2] is linearized and split in the form w' = {A + B)w , the solution over a timestep 
T requires an approximation to the operator e T ( A + B ) _ High-order approximations may be obtained 
through appropriate choice of partial steps Tj where 


= e k J e kj - 1 ' 


(25) 


Formally, with support from numerical studies 37, 5l|, the splitting scheme given by Eq. 1251 may 
be may be extended to the semi-linear parabolic form of Eq. [2] given by 


w 


= Aw + fs(w) 


(26) 


by replacing the exponential operator e Tk i B with a step of the nonlinear equation w' = /b (w) over 
the interval T\ . For reference, the complex splitting schemes used in this work are provided in 
Table HOI 


4-2. Brusselator 

The Brusselator [52, [2l| is a stiff nonlinear diffusion-reaction problem describing chemical kinet¬ 
ics of a tri-molecular chemical reaction. The test case considered here is a two-dimensional hybrid 
of the one- and two-dimensional Brusselator problems presented by Hairer et al. [HJ, and Hairer 
and Wanner [3l| , with governing equations given by 


dv/dt = e (d 2 v/dx\ 2 + d 2 v/dx 2 2 ) + A — (B + l)v + wv 2 , 

dw/dt = e [d 2 w/dx\ 2 + d 2 w/dx 2 2 )+Bv — v 2 w, (27) 

and initial conditions u(0, x ) = A+sin(27ra:), i>(0, x) = B/A+cos(2iry). The initial state is therefore 
a simple perturbation of the equilibrium solution. The problem is configured with parameters 

e = 0.02, A = 1, and B = 3, and the solution is obtained at t = 2, or t = 8, on the domain 

0 <a:i<l, 0 <X 2 <l, under periodic boundary conditions. 
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Li error (species v) 

Figure 4: Li errors plotted against number of timesteps, Nt, for species v of the two-dimensional 
Brusselator problem. Results correspond to the split problem, with linear diffusion treated via frkc 
methods at orders 2, 4, and 6, and nonlinear reaction terms integrated via standard techniques. 
Guide lines are shown for (Lierror) -1 / 2 , (Lierror) -1 / 4 , (Lierror)” 1 / 6 . Table [T] gives the values for 
all points shown. 
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Table 1: Error convergence results for the two-dimensional Brusselator test problem solved via split 
frkCs schemes with TV = 2, 4, 6. Each row corresponds to a specific test with columns listing: TV, 
order of accuracy; TVV, the number of timesteps; Li norm of the error between the approximate 
and exact solutions; Li order of convergence with reference to previous row; Loo error; Loo order. 
Errors refer to the solution for species v. Li errors are also shown in Fig HI 


TV 

Nt 

Li error 

Li order 

Loo error 

Loo order 

2 

50 

1.92 x 10~ 4 

- 

4.25 x 10~ 4 

- 


100 

4.76 x 10~ 5 

2.02 

1.03 x 10~ 4 

2.04 


200 

1.18 x 10~ 5 

2.01 

2.55 x 10~ 5 

2.02 


400 

2.95 x 10~ 6 

2.00 

6.32 x 10~ 6 

2.01 


800 

7.36 x 10~ 7 

2.00 

1.57 x 10~ 6 

2.01 


1600 

1.83 x 10- 7 

2.01 

3.92 x 10- 7 

2.01 


3200 

4.53 x 10“ 8 

2.02 

9.68 x 10“ 8 

2.02 


6400 

1.08 x 10“ 8 

2.07 

2.30 x 10“ 8 

2.07 


12 800 

2.16 x 10~ 9 

2.32 

4.61 x 10~ 9 

2.32 

4 

50 

7.85 x 10~ 6 

- 

1.21 x lO^ 5 

- 


100 

1.71 x 10- 7 

5.52 

2.76 x 10- 7 

5.46 


200 

9.96 x 10~ 9 

4.10 

1.65 x 10“ 8 

4.06 


400 

6.05 x 10- 10 

4.04 

1.02 x lO^ 9 

4.02 


800 

3.76 x 10- 11 

4.01 

6.36 x 10~ n 

4.00 


1600 

2.35 x 10~ 12 

4.00 

3.98 x 10 -12 

4.00 


3200 

1.47 x 10- 13 

4.00 

2.49 x 10~ 13 

4.00 


6400 

9.13 x 10~ 15 

4.01 

1.55 x lO^ 14 

4.00 


12 800 

5.37 x 10- 16 

4.09 

9.15 x KT 16 

4.09 

6 

50 

4.41 x 10~ 5 

- 

7.37 x 10~ 5 

- 


100 

9.92 x 10- 7 

5.48 

2.35 x lO^ 6 

4.97 


200 

1.31 x 10“ 8 

6.24 

2.88 x 10“ 8 

6.35 


400 

1.04 x 10- 10 

6.98 

1.78 x lO^ 10 

7.34 


800 

1.39 x 10~ 12 

6.23 

2.34 x 10 -12 

6.25 


1600 

2.06 x 10~ 14 

6.07 

3.52 x lO^ 14 

6.05 


3200 

3.16 x 10~ 16 

6.02 

5.49 x lO^ 16 

6.01 


6400 

4.44 x 10~ 18 

6.16 

9.87 x KT 18 

5.80 


12 800 

7.01 x 10~ 19 

2.66 

2.01 x lO^ 18 

2.30 
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4-3. Linear order conditions 

The semi-discrete form of Eq. [57] may be written w' = Aw + /s(tc), where A describes the 
discretization of the Laplacian with respect to x\ and * 2 , and /b(w) contains the reaction terms. 
Linear diffusion terms are integrated using frkc methods and nonlinear reaction terms via stan¬ 
dard techniques. The linear order properties of the frkc stability polynomials are confirmed by 
considering the convergence rates of the approximated solution to the exact solution at t = 2 as a 
function of step size. 

For all presented results, we use M = 20, and the approximation 7 m = 1- The number of 
grid points is 400 in each of the two spatial variables. For these parameters, the frkc stability 
polynomials achieve approximately 81% (/ 3r = 1066.667), 74% (/3# = 1600), and 73% (/ 3 _r = 
2133.333) of the optimal intervals for TV = 2, 4, 6 respectively. With respect to the corresponding 
standard explicit Runge-Kutta schemes, these values represent a speedup in efficiency by factors of 
approximately 27 for TV = 2, and 30 for both TV = 4 and TV = 6 . All polynomials are damped with 
damping parameter vq = 0.05, reducing the stability domains’ real extents by factors of 1 — z'o/TV. 
Finally, in order to meet the specified solution time, timesteps are scaled by 0.9846, 0.9001, 0.7563 
for TV = 2, 4, 6 respectively. Quadruple precision is used in all calculations. Results are presented in 
Table [T] where the Li and Lqo errors are shown over a range of resolutions at each considered value 
of TV. Fig. [4] illustrates the dependence of the Li errors on the number of timesteps, TVt, for species 
v. With the exception of the final point for the sixth-order integration, where machine precision 
is exceeded, all solutions are converging in good agreement with the nominal orders of accuracy 
(i.e. (error) _1 / Ar ). Fitting the Li errors yields observed orders 2.04 ± 0.01, 4.08 ± 0.04, 6.1 ± 0.2 
for TV = 2, 4, 6 respectively, while the errors give 2.05 ± 0.01, 4.09 ± 0.05, and 6.0 ± 0.2. We 
conclude that frkc methods demonstrate internal stability and comply with linear order conditions 
to the specified order of accuracy. 

4-4- Second-order comparative studies 

Since all order conditions are linear at second-order, FRKC2 schemes will naturally maintain 
second-order accuracy for nonlinear problems without the necessity of splitting or composition 
methods. Here we present comparative studies between FRKC2 and a number of alternative numer¬ 
ical integration methods. In particular, we provide comparisons with the rkc method 0 which, 
similarly to FRKC2, depends on the properties of Chebyshev polynomials. We also compare results 
with a GMRES Krylov-preconditioned BDF integrator from the cvode numerical integration pack¬ 
age [53|. The cvode solver maintains a specified tolerance by means of adaptive stepping up to a 
maximum fifth-order accuracy. However, the order is restricted to 2 for the CVODE2 solver used in 
these comparisons. 

We proceed by considering the two-dimensional Brusselator problem described in Section 14.21 
with the solution taken at time t = 8. The stepsize is fixed for individual tests of the explicit 
schemes and the number of internal stages is optimized for the selected stepsize. As such, each of 
the numerical solutions generated for these tests is derived from a single distinct stability polyno¬ 
mial. In general, however, error control procedures may be implemented 0 which will result in 
stability polynomials of varied degree contributing to particular solutions. The optimal efficiency 
for extended stability explicit solvers follows Iwall c* (error) _1 / 2Ar (where Xwall is the wall-time 
required for computation of a particular solution). 

Results are provided in Table [2] for FRKC2, rkc, and CVODE2. The Li errors for species v are 
plotted in Fig. [5] (a) against the time required for the simulations to be carried out on a standard 
desktop machine at double precision. While the FRKC2 solver requires complex arithmetic, this is 
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Table 2: Errors from FRKC2, rkc, CVODE2 from tests of the two-dimensional Brusselator problem. 
The number of timesteps, Nt, and the number of stages per timestep, L , (or the error tolerance, 
Tol, in the case of cvode) are given in the first two columns respectively. The wall time taken for 
each run, 7w all, is presented in the third column. Li and Loo errors for both species are presented 
in the remaining columns. The Li error for species v is plotted in Fig. [5] (a). 


N t 

L /Tol 

Twall (s) 

Species v 

Li error Loo error 

Species w 

Li error L^ error 

FRKC2 

50 

80 

25 

5.16 x 

10~ 3 

5.52 x 

10~ 3 

1.30 x 

10~ 3 

1.40 x 10~ 3 

137 

48 

42 

1.41 x 

10~ 3 

1.45 x 

10~ 3 

6.06 x 

10~ 4 

6.25 x 10~ 4 

308 

32 

63 

3.24 x 

10~ 4 

3.34 x 

10~ 4 

1.64 x 

10~ 4 

1.68 x 10~ 4 

548 

24 

84 

1.08 x 

10~ 4 

1.11 X 

10~ 4 

5.66 x 

10~ 5 

5.81 x 10~ 5 

1317 

16 

137 

1.85 x 

10~ 5 

1.90 x 

10~ 5 

1.00 x 

10~ 5 

1.03 x 10~ 5 

2341 

12 

186 

6.03 x 

10~ 6 

6.20 x 

10~ 6 

3.30 x 

10~ 6 

3.38 x 10~ 6 

5266 

8 

288 

1.26 x 

10~ 6 

1.30 x 

10~ 6 

6.94 x 

io- 7 

7.11 x 10- 7 

9361 

6 

399 

4.23 x 

io- 7 

4.35 x 

io- 7 

2.32 x 

io- 7 

2.38 x 10- 7 

21062 

4 

637 

9.48 x 

10“ 8 

9.76 x 

10~ 8 

5.16 x 

10~ 8 

5.30 x 10 8 

05 026 

2 

1881 

5.71 x 

10~ 9 

5.88 x 

10~ 9 

3.11 x 

10~ 9 

3.19 x 10~ 9 

RKC 

39 

90 

32 

4.97 x 

10~ 3 

5.81 x 

10~ 3 

8.05 x 

10~ 3 

8.54 x 10~ 3 

78 

64 

45 

6.59 x 

10~ 4 

7.64 x 

10~ 4 

6.50 x 

10~ 4 

7.02 x 10~ 4 

137 

48 

60 

2.49 x 

10~ 4 

2.83 x 

10~ 4 

1.08 x 

10~ 4 

1.26 x 10~ 4 

309 

32 

90 

5.07 x 

10~ 5 

5.73 x 

10~ 5 

1.36 x 

10~ 5 

1.75 x 10~ 5 

549 

24 

119 

1.59 x 

10~ 5 

1.81 x 

10~ 5 

3.84 x 

10~ 6 

5.12 x 10~ 6 

1237 

16 

177 

3.13 x 

10~ 6 

3.55 x 

io~ 6 

7.27 x 

10- 7 

9.85 x 10~ 7 

2206 

12 

237 

9.92 x 

10- 7 

1.13 x 

10~ 6 

2.28 x 

10- 7 

3.10 x 10- 7 

5007 

8 

359 

2.00 x 

10- 7 

2.27 x 

10- 7 

4.50 x 

10~ 8 

6.17 x 10~ 8 

9012 

6 

486 

6.57 x 

10“ 8 

7.46 x 

10 8 

1.40 x 

10~ 8 

1.94 x 10~ 8 

21028 

4 

750 

1.52 x 

10“ 8 

1.71 x 

10~ 8 

1.88 x 

10~ 9 

3.03 x 10~ 9 

39 428 

3 

1056 

6.56 x 

10~ 9 

7.20 x 

10~ 9 

4.13 x 

IO- 10 

7.98 x 10~ 10 

CVODE2 

1226 

5 x 10~ 6 

107 

2.37 x 

10~ 3 

2.42 x 

10~ 3 

9.65 x 

10~ 4 

1.01 x 10~ 3 

1499 

10~ 6 

128 

7.79 x 

10~ 4 

8.26 x 

10~ 4 

5.85 x 

10~ 5 

8.45 x 10~ 5 

2534 

io- 7 

194 

4.69 x 

10~ 5 

5.83 x 

10~ 5 

2.99 x 

10~ 5 

3.91 x 10~ 5 

4991 

10~ 8 

312 

3.06 x 

10~ 5 

3.22 x 

10~ 5 

3.13 x 

10~ 5 

3.27 x 10~ 5 

10 029 

10~ 9 

523 

5.47 x 

10~ 6 

5.60 x 

10~ 6 

4.35 x 

10~ 6 

4.45 x 10~ 6 

22 763 

10- 10 

987 

7.05 x 

10- 7 

7.14 x 

10- 7 

5.06 x 

10- 7 

5.08 x 10~ 7 

48 444 

10- 11 

1706 

1.55 x 

10- 7 

1.58 x 

10- 7 

1.07 x 

10- 7 

1.09 x 10- 7 

09474 

10- 12 

3405 

2.96 x 

10“ 8 

3.01 x 

10“ 8 

2.06 x 

10~ 8 

2.10 x 10~ 8 

132430 

10- 13 

6756 

6.05 x 

10~ 9 

6.16 x 

10~ 9 

4.16 x 

10~ 9 

4.24 x 10~ 9 
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Figure 5: Performance results derived from the Li error in the solution to species v in the stiff 
nonlinear Brusselator problem. Panel (a) shows results for the second-order schemes FRKC2, rkc, 
and CVODE2 (see Table A guide line is shown for (Lierror) 1 / 4 . Panel (b) shows data for the 
higher-order schemes frkC4s, frkC6s, ROCK4, and cvode (see Table m. Guide lines are shown for 
(Lierror) -1 / 8 and (Lierror) -1 / 12 . 


compensated by smaller errors than for the rkc solver at equivalent numbers of timesteps. Overall, 
FRKC2 runs at about 70% of the efficiency of rkc. As previously noted, following a similar strategy 


to Lebedev [111 . I28j will improve performance. 


4-5. High-order comparative studies 

An advantage of frkc methods over other extended stability methods is extensibility to arbi¬ 
trarily high-order linear stability polynomials. In order to apply high- (above second-) order frkc 
stability polynomials to nonlinear problems, complex splitting techniques may be employed (as 
demonstrated in Sec. 14.311 . In the following tests, we consider fourth- and sixth-order solutions 
of the two-dimensional Brusselator problem. We note that finishing stages based on the theory 
of the composition of Butcher series may also be may be used to meet nonlinear order condi¬ 
tions |2ll Il2t 17]. While composition methods may, in principle, offer improved efficiency over 
splitting techniques, in the case of ROCK4, the application of finishing stages has been observed 
to result in order reduction problems, as well as erratic convergence prop erties, limiting the num¬ 
ber of internal stages to a relatively small number of internal stages [2J, |35[. (The limit adopted 
within the ROCK4 code is L = 152.) Furthermore, the number of nonlinear order conditions, and 
hence the complexity of the composition strategy, grows rapidly with increasing order [22|: there 
are four nonlinear order conditions at fourth-order, 31 order conditions at sixth-order, and 192 at 
eighth-order. 
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Table 3: Errors from frkC4s, frkC6s, ROCK4, cvode from tests of the two-dimensional Brusselator 
problem. The number of timesteps, Nt, and the number of stages per timestep, L , (or the error 
tolerance, Tol, in the case of cvode) are given in the first two columns respectively. The wall time 
taken for each run, Twall, is presented in the third column. Li and Loo errors for both species are 
presented in the remaining columns. The Li error for species v is plotted in Fig. El (b). 





Species v 

Species w 

TVj 1 

L /Tol 

'■L WALL (s) 

Li error 

Lqo error 

Li error 

Lqo error 

FRKC4s 

23 

96 

31 

8.27 x 10~ 4 

8.49 x 10~ 4 

5.49 x 10~ 4 

5.62 x 10~ 4 

51 

64 

48 

3.30 x 10~ 5 

3.30 x 10~ 5 

2.35 x 10~ 5 

2.37 x 10~ 5 

91 

48 

70 

3.21 x 10~ 6 

3.25 x 10~ 6 

2.31 x 10~ 6 

2.32 x 10~ 6 

222 

32 

117 

9.17 x 10“ 8 

9.27 x 10~ 8 

6.66 x 10“ 8 

6.68 x 10~ 8 

395 

24 

166 

9.17 x 10~ 9 

9.27 x 10~ 9 

6.68 x 10~ 9 

6.71 x 10~ 9 

887 

16 

277 

3.34 x 10~ 10 

3.38 x 10- 10 

2.42 x 10- 10 

2.43 x 10- 10 

1577 

12 

411 

7.83 x 10~ 12 

7.91 x IQ- 12 

4.03 x 10~ 12 

4.09 x IQ- 12 

FRKC6s 

5 

144 

41 

1.00 x 10~ 3 

1.02 x 10~ 3 

8.60 x 10~ 4 

8.74 x 10~ 4 

10 

96 

56 

4.58 x 10~ 5 

4.62 x 10~ 5 

3.57 x 10~ 5 

3.61 x 10~ 5 

17 

72 

75 

2.26 x 10~ 6 

2.31 x 10~ 6 

1.79 x 10~ 6 

1.84 x 10~ 6 

42 

48 

125 

1.97 x 10~ 8 

2.01 x 10~ 8 

1.55 x 10“ 8 

1.57 x 10~ 8 

75 

36 

184 

1.04 x 10~ 10 

1.04 x 10- 10 

8.27 x 10- 11 

8.33 x 10- 11 

168 

24 

296 

1.70 x lO^ 11 

1.72 x 10~ n 

1.34 x IQ- 11 

1.35 x 10~ n 

ROCK4 

56 

102 

16 

1.31 x 10~ 3 

7.26 x 10~ 3 

1.34 x 10~ 3 

7.19 x 10~ 3 

130 

67 

25 

3.79 x 10~ 5 

2.13 x 10~ 4 

2.88 x 10~ 5 

1.49 x 10~ 4 

224 

51 

32 

3.40 x 10~ 6 

1.83 x 10~ 5 

2.76 x 10~ 6 

1.42 x 10~ 5 

451 

36 

46 

1.64 x 10- 7 

8.18 x 10- 7 

1.38 x 10- 7 

7.56 x 10- 7 

749 

28 

59 

1.95 x 10“ 8 

8.90 x 10~ 8 

1.60 x 10“ 8 

8.14 x 10~ 8 

1483 

20 

84 

1.14 x 10~ 9 

4.64 x 10~ 9 

8.99 x 10- 10 

4.74 x 10~ 9 

2345 

16 

107 

1.65 x 10- 10 

5.58 x 10- 10 

1.25 x 10- 10 

5.66 x 10- 10 

4282 

12 

146 

6.48 x 10~ 12 

3.32 x 10~ n 

6.13 x 10~ 12 

2.80 x 10~ n 

CVODE 

2001 

10~ 6 

151 

1.48 x 10~ 3 

1.51 x 10~ 3 

1.13 x 10~ 3 

1.15 x 10~ 3 

3203 

io- 7 

209 

4.38 x 10- 7 

1.26 x 10~ 6 

4.83 x 10- 7 

1.71 x 10~ 6 

3309 

10~ 8 

250 

3.74 x 10~ 5 

3.90 x 10~ 5 

2.91 x 10~ 5 

3.03 x 10~ 5 

6628 

10~ 9 

407 

1.09 x 10~ 9 

5.14 x 10~ 9 

9.33 x 10- 10 

4.06 x 10~ 9 

4818 

10-io 

307 

3.88 x 10- 7 

3.89 x 10~ 7 

3.27 x 10- 7 

3.28 x 10- 7 

6246 

io- 11 

416 

8.46 x 10~ 9 

8.97 x 10~ 9 

4.88 x 10~ 9 

4.95 x 10~ 9 

8376 

io- 12 

503 

1.69 x 10~ 9 

1.76 x 10~ 9 

1.32 x 10~ 9 

1.37 x 10~ 9 

10515 

io- 13 

584 

1.96 x 10- 10 

1.97 x 10- 10 

1.38 x 10- 10 

1.41 x IO- 10 

16 036 

10 14 

817 

4.79 x lO^ 12 

5.29 x 10~ 12 

3.38 x 10~ 12 

3.84 x IQ- 12 


17 












We present comparisons of the split schemes frkC4s and frkC6s with the fourth-order ROCK4 
scheme. Reference solutions obtained using the cvode solver are also presented for integrations 
carried out to a maximum fifth-order accuracy. The test conditions are otherwise as described 
in Sec. 14.41 All data are presented in Tabic El with Li errors for species v plotted in Fig. [5] (b) 
against the time required for the simulations at double precision. We note that frkC4s is shown 
to run at approximately half the efficiency of ROCK4. This is primarily due to the additional 
computational overhead of carrying out calculations with complex values quantities, which, as 
noted, may be countered to some degree by rolling conjugate pair calculations together using the 
scheme presented by Lebedev (TT], [28[ . In terms of efficiency, for the presented problem, the frkC6s 
and frkC4s methods lie approximately midway between ROCK4 and cvode. Except at the very lowest 
acceleration parameters considered, the frkCs trials show the predicted behavior (i.e. Iwall oc 
(Lierror) -1 / 2Ar ). 

5. Conclusions 

The fully prescribed analytic form of a new class of extended stability polynomials which sat¬ 
isfy all required linear order conditions to arbitrarily high-order has been presented. Factorized 
Runge-Kutta-Chebyshev (frkc) stability polynomials are derived from first principles by inductive 
considerations of the implied recurrence relations. At order N, the frkc polynomial of rank N, 
and degree L = MN , is shown to have the form of a summation of Chebyshev polynomials, with 
degrees at intervals of M, up to degree L. The N + 1 weightings of the contributing Chebyshev 
polynomials are chosen to comply with the N linear order conditions, coupled with a constraint. 
A damping procedure for broadening the stability domain of the frkc stability polynomials to a 
finite width along the real axis is described which preserves the order of accuracy. The resultant 
stability polynomials have been demonstrated to have 81%, 74% and 73% of the optimal intervals 
for orders 2, 4, 6 respectively, frkc numerical integration schemes are represented as a sequence 
of L sequenced forward Euler steps (stages) involving complex-valued timesteps constructed from 
the roots of frkc stability polynomials of degree L. Internal stability is maintained by means of 
a sequencing algorithm, which limits the maximum internal amplification factor to ~ L 2 : reserv¬ 
ing 8 digits for accuracy, a hypothetical scheme of 10,000 stages is therefore viable in a numerical 
integration carried out at 16 digit precision. 

Split frkCs schemes have been applied at orders 2, 4, and 6, to the linear diffusion operator 
in numerical experiments on a stiff two-dimensional Brusselator reaction-diffusion system leading 
to the verification of expected convergence rates, and hence compliance with the necessary linear 
order conditions. 

We have presented comparative studies of the performance of FRKC2 with rkc, an established 
explicit extended stability code, and CVODE2, an implicit preconditioned BDF solver from the cvode 
suite limited to second-order accuracy. FRKC2 has been shown to be substantially more efficient 
than the CVODE2 solver, while performing at about 70% of the efficiency of rkc. 

At higher orders, nonlinear order conditions require special attention. We have considered 
treatment of these nonlinear conditions through complex splitting techniques in efficiency tests 
of higher order (4 and 6) split frkCs schemes in comparison with results from the the fourth- 
order ROCK4 code, which uses composition methods, and the implicit fifth-order cvode solver. 
The tested frkCs methods are found to have intermediate efficiency to ROCK4 and cvode . We 
propose implementing conjugate pairing and Butcher series composition methods in future high- 
order implementations of frkc methods. 
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Appendix A. Scheme patterns M = 20 


d 


2 _ 
0 — 


d 


2 _ 

1 


(I 


2 _ 
2 — 


267 

400 

1 

1800 

1201 

7200 


(A.l) 


3126039467 

6144000000 

244573733 

7680000000 

3212226667 

15360000000 

63194381 

7680000000 

789861181 

61440000000 


(A.2) 


7446093942631413209 

17915904000000000000 

158532158867283313 

2985984000000000000 

1022936325403301087 

4777574400000000000 

35821864811075087 

10749542400000000000 

1048968349471238687 

35831808000000000000 

32100268736824717 

17915904000000000000 

180240686854539517 

214990848000000000000 


(A.3) 


23 
















Appendix B. Splitting schemes 


Table B.4: Complex operator splitting parameters for N = 2, 4, 6 [HU, H3, HU. The final row for 
each quoted value of N lists: J, the number of distinct sweep configurations required; ki ■ ■ • kj , 
the sequence of J sweeps, labeled by j , required for a single extended interval to order-JV. The 
remaining rows are in pairs listing: j, the index of the distinct sweep; s ft(2j), the real component 
of the sweep timescale; 3(Tj), the second row lists the imaginary part of the sweep timescale. 


N 

J 

HCZi) 

3W) 


J 

fci • • • kj 

2 

1 

1.0 

0.0 


2 

0.5 

0.0 


3 

212 

4 

1 

1/4 

n 


2 

1/10 

-1/30 


3 

4/15 

2/15 


4 

4/15 

-1/5 


9 

213141312 

6 

1 

0.0625 

0.0 


2 

0.024 694 876 087 018 064 640 910 864 996 842 247 838 60 
-0.007 874 795 562 906 877 058171577 949 526 942163 20 


3 

0.063 813 474 021302 699 779 366 304188 200146 963 20 
0.035 365 761034143 327 804 629 404 649 714 741812 70 


4 

0.068 425 094 030 316 441970 397 007 821744 684 058 50 
-0.062 262 244 450 748 676 995 332 540 644 447 596 04610 


5 

0.088 047 701092 267 837 626 997195 869 408 667 577 20 
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